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ABSTRACT 

During the formation of the Milky Way, > 100 central black holes (BHs) may have 
been ejected from their small host galaxies as a result of asymmetric gravitational wave 
emission. We previously showed that many of these BHs are surrounded by a compact 
cluster of stars that remained bound to the BH during the ejection process. In this 
paper, we perform long term A'^-body simulations of these star clusters to determine the 
distribution of stars in these clusters today. These numerical simulations, reconciled 
with our Fokker-Planck simulations, show that stellar density profile follows a power- 
law with slope ~ —2.15, and show that large angle scattering and tidal disruptions 
remove 20 — 90% of the stars by ^ 10^° yr. We then analyze the photometric and 
spectroscopic properties of recoiled clusters accounting for the small number of stars 
in the clusters. We use our results to perform a systematic search for candidates 
in the Sloan Digital Sky Survey. We find no spectroscopic candidates, in agreement 
with our expectations from the completeness of the survey. Using generic photometric 
models of present day clusters we identify ^ 100 recoiling cluster candidates. Follow-up 
spectroscopy would be able to determine the nature of these candidates. 

Key words: galaxies:kinematics and dynamics-galaxies: nuclei-black hole physics- 
gravitational waves-star clusters 



1 INTRODUCTION 
1.1 Motivation 

In the standard cosmological context of hierarchical galaxy 
formation, the first galaxies to form were small (-^ 10* 
and g rew through subsequent accretion and mergers (|Loebl 
I2OIOI ). If central black holes (BHs) were common in the 
earliest epochs of galaxy formation, then for most major 
mergers, the two BHs would also eventually merge. If there 
were any asymmetry in the inspiral and coalescence of two 
BHs, whether due to a difference in mass, or the alignment 
of the BHs' spin vectors, then there would inevitabl y be a 
net l i near momentum k ick to the mer ger remnant l|Peresl 
1 19621 : iBekensteinI 1 19731 : iFitchettI ll983D . This kick, which 
is typically h undreds of kms~ ^ for mergers with comp a- 
rable ma sses ("Baker et al. l l2006l : ICampaneni et al. I l2007bl lal: 
iTichv fc' M arronctti 20o3), is independent of the total mass 



merger in the first galaxies would inevitably lead to the 
ejection of the BH. Interestingly, the typical kick velocity 
is smaller than the escape velocity of the Milky Way halo, 
and so, although the first galaxies to merge would have 
lost their BHs, many of these 'r ogue' BHs should remain 
in the Milky Way halo today (iMadau fc QuataertI |2004|: 
Volonteri fc Pernall2005l: iLibeskind et al.ll2006l : iMicic et all 
20081 : lO'Leary fc Loebll2009l v 



Before coalescence, the BH-BH binary may be sur- 
rounded by both a disk of dense gas and a dense cusp 
of stars. For bound matter with orbital velocities larger 
than the kick velocity, the gravitational wave kick per- 
turbs the orbit of the material, but does not unbind it 
from the recoiled BH even if the BH is ejected from the 
galaxy fLocb 2007). For gas disks, viscosity eventually 
causes the BH to accrete the surrounding gas on or der of 
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cal effect in the smallest galaxies (iMadau & Ouataertj 
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Schnittman 20071: Blecha & Loebl 


20081). Indeed, a maior 
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Myr, leavmg me ati as a snort lived quasar (IL/oet 
Blecha fc Loebll2008l : iGuedes et allboiol : ISiiacki et al' 



Blecha et al.ll201ll ). After depleting all bound gas, the 



BHs will only be visible if they pas s through dense g as 
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in the galaxy and reacc r ete material (Islam et al.l 12004 



Volonteri fc Pernal [20051: iMii fc Totanirf2005l : iBertone et all 
2005; Blec ha fc Loebll2008l ). Stars, on the other hand, are ef- 
fectively collisionless, and will remain bound as a long lived 
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system JCualandris fc MerrTt3 20081: Komossa fc Merrit3 



I2OO8I : lO'Learv fc Loebll2009l : iMerritt et al.ll2009l ). and may 
actual ly be tidally disr upted by the BH after ejec- 
tion (iKomossa fc MerrittI 12008) or produce winds pro- 
viding a new source of gas accretion to the black hole 
ijO'Learv fc Loebll2009l l. 

In O'Leary fc Loeb (2009; hereafter as Paper I), we used 
> 1000 merger tree histories of the Milky Way galaxy to cal- 
culate the number and mass distribution of these recoiled 
BHs. We assumed that after a major merger, with galaxy 
mass ratio greater than 1:10, the black holes merged imme- 
diately. We found that > 100 BHs with M. > lO"* M© should 
be in the halo today, surrounded by c ompact star clusters 
that are ~ 1% of their BH's mass fsee lMerritt et alJlioogl . 
for a similar discussion applied to the nearby Virgo Cluster). 
These clusters are expected to roughly follow the dark mat- 
ter distribution of the halo, since they have kick velocities 
typically less than the velocity dispersion of the galaxy. At 
their typical distances dynamical friction is not important 
over a Hubble time. The most massive star clusters have a 
much higher internal velocity dispersion than globular clus- 
ters because they are gravitationally bound by the BH at 
their center. Finding these clusters in the Milky Way halo 
will allow us to peer back in time and look at some of the first 
BHs to form. In this paper we re-investigate the long term 
evolution of recoiled star clusters using full A''-body simula- 
tions, with a one-to-one correspondence between stars and 
A'^-body particles. We also extend the Fokker-Planck simula- 
tions from Paper I to include large-angle scattering between 
stars in order to reconcile these new simulations with our 
results from Paper I, and extend the reach of our simula- 
tions to larger BHs. We use the results of these simulations 
to generate the photometric properties of recoiled clusters. 
Because these clusters have so few stars, stochastic variation 
dominates over the dispersion of the cluster colors. Instead 
of using averaged stellar evolution tracks of old star systems, 
we use a Monte-Carlo approach to generate individual star 
cluster colors and sizes to identify the typical properties of 
these systems and compare them to the properties of stars 
and galaxies in the the Sloan Digital Sky Surve^ D ata Re- 
lease 7 (hereafter SDSS DRT lAbazaiian et aLlfjoogj ). 

Our paper is organized as follows. In § 11.21 we describe 
recoiled clusters and briefly summarize the main results from 
Paper I. In § [2l we use A-body simulations to follow the 
long term dynamical evolution of recoiled star clusters. We 
then extend our previous Fokker-Planck analysis in § [3] to 
include the ejection of stars from large angle scattering. In 
§ m we develop a series of simple photometric models that 
we use to search for recoiled clusters in §(5] In §[6l we search 
the spectroscopic database of SDSS for candidate clusters. 
Finally in § [HI we summarize our paper and describe the 
main conclusions. 



1.2 Stellar Cusps and Recoiled Black Holes 

For a relaxed stellar s ystem around a central massive BH, 
iBahcall fc Wolj (|l976l ') found that the stellar density follows 
a power-law profile within the radius of influence of the BH, 
Ti — GM,/al, where M, is the BH mass, and cr* is the 



stellar velocity dispersion. Out to the radius where the total 
mass in stars around the BH is twice the mass of the BH 
the density profile is 



M. 3- 



2'Krf 



(1) 



where a = 1.75 for a cluster of single mass stars of mass m*. 
If the binary black hole merge s on a timescale compa rable 
to the relaxation time, then the lBahcall fc Wolj (|l976l ) cusp 
will be regenerated as the binary inspirals. However, if the 
binary merges on a shorter timescale than the relaxation 
time (e.g., through its interaction with gas ), then the stel- 
lar density profile is expected to be much shallower with 
Q ~ 1 and with few er stars within (jMerritt fc Sze Ul l2006l : 
iMerritt et"aLl 120071 '). For the BH masses we consider here, 
the relaxation timescale is much less than the age of the 
universe, 



tr{r) 



10" yr 
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(2) 



The kick on the BH remnant occurs over a timescale 
much shorter than the orbital time of the stars. Therefore, 
in the frame of the BH, the stars all instantaneously receive a 
kick, — Hk. To first order, stars with a total energy < — m^w^, 
will remain bound to the BH as it is ejected from the galaxy. 
For the Keplerian potential of the BH, this approximately 
corresponds to all stars with r < rk = (wk/o"*)~^ri. From 
equation ([T}, this corresponds to a immber of stars 



Ad 



2M. 



Vk 



: 4x10'' 



M. 



105 M0 



Wk 
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where we set a = 1.75, and used the M, 



5.6cr* 

(3) 

a-i, relation to 



^ |http : //www ■ sdss ■ org/dr7/ | 



determine ri (|Tremaine et al.ll2002f ). 

The star cluster will begin to expand away from the BH 
immediately after it is ejected from the galaxy. In Paper I, we 
followed the long term evolution of the star cluster by nu- 
merically integrating the time dependent, one-dimensional 
Fokker-Planck equation for stars around a central BH. We 
found that the density cusp of stars around the BH quickly 
expands to the point that the relaxation timescale of the 
system U is approximately its age tn- In our simulations, 
the total mass in stars was roughly constant. However, our 
calculations could not consistently account for either strong 
encounters or resonant interactions between stars. In both 
cases, we would expect a larger fraction of stars to be ejected 
or destroyed. 

We estimated the number of recoiled BHs in the Milky 
Way Halo using an ensemble of > 10^ merger tree histo- 
ries of the Milky Way, con volved with the probability dis - 
tribution of kick velocities (|Schnittman fc Buonannoll2007l '). 
We found that a typical Milky Way like galaxy contains as 
many as 100 recoiled BHs with M, > IO^Mq, with a mass 
spectrum dA/dM. (x M^^ . Because the kick velocity distri- 
bution peaks at low velocities, ~ lO^kms"'^, the majority 
of recoiled clusters had the minimal kick velocity needed to 
escape from the host galaxy. In Paper I, we estimated this 
to be Vbsc ~ 5.6(7* immediately after the galaxy merger. 
Overall, these results were consistent with previous stud- 
ies that looked at a population of BHs in the Milky Way 
halo, whether from gravitational wave recoil, through three 
body encounters, or as the remnants of the seed population 
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of black holes (|Madau fc Quataert 2004 : Volpnteri fc Pernal 



I2OO5I : iLibeskind et al.ll2006l : iMicic etal 



2008h 



2 iV-BODY SIMULATIONS 

Small star clusters around recoiled BHs present an interest- 
ing dynamical system that can be modeled directly in A^- 
body simulations, as well as through approximate methods 
such as solving the Fokker-Planck equation. Because some 
star clusters have only a few thousand stars, it is possible to 
simulate the star clusters with a one to one correspondence 
between stars and A'^-body particles. In this section, we sim- 
ulate the star clusters directly, and compare the results to 
our previous Fokker-Planck simulations from Paper I. 



2.1 Method 

We use the publicly availab le A'^-body code BHintQ 
iLockmann fc BaumgardU |2008|) to simulate the long term 
evolution of star clusters around recoiled BHs. BHint was 
developed to precisely integrate the equations of motion of 
stellar systems around massive BHs, where the BH domi- 
nates the motion of the stars. 

The initial conditions for the re coiled BH and st a r clus - 
ter are set up by first generating iBahcall fc WoUl (Il976l ) 
stellar density cusps around a M, = 10* BH, follow- 



ing equation ((T}, with n* 



such that the total mass 
of stars within Vi is twice that of the black hole. The stars' 
velocities are initially selected from a Gaussian distribution 
with crj — GM,/r / {I + a), where a — 1.75, which is an 
excellent approximation to the velocity distribution of stars 
in a power-law density cusp. We use two model mass func- 
tions for the stars. In Model I we use equal-mass stars with 
JTiv, = 1 Mq . In Model II, we use a more realistic mass func- 
tion to model a population of old stars. Although massive 
stars likely play an important role in the evolution of clus- 
ters, they are short lived compared to the cluster lifetime. 
Stars with mass < 20 M©, evolve to form neutron stars or 
white dwarfs, which are comparable in mass to long lived 
main sequence stars, and so should not significantly alter 
the dynamics of the system except to increase the mass- 
to-light ratic0. More massive stars form black holes, with 
mass ~ 10 M©, which may dramatically alter the dynamics 
of the system. The se black holes may even segregate before 
the binary merges l|Morris| [l993: Mcrritt 2009). If the black 
holes dominate the dynamical evolution of the system, such 
systems may not have any luminous stars to observe. How- 
ever the fraction of black holes in the region immediately 
after reforming the cusp is highly uncertain, and so we take 
the extreme opposite approach and ignore the black hole 
population. Assuming that there are very few black holes, 
which may be the case in a subset of the recoiled clusters, 
we use a relatively flat mass function for low mass stars 
(diV/dm* oc m~^-^^) with 0.2 Mq < m* < 1.0 M©. The 
mass of each star is generated randomly following this dis- 
tribution until the total mass of the cusp reaches 2M,. We 



^ Available at http : / /www . astro . uni-bonn . de/english/downloads 

^ The total amount of mass lost is < < A/, , and will accordingly 
only perturb tlie orbits of the stars. 



use such a shallow power law bec ause of the b reak in the 
initial mass function at « 0.5 Mq (|Kroupall200"ll V 

After generating the cusps for Models 1 fc 2, we then 
kick each star with a velocity Vk — VkZ, and remove all stars 
that are unbound to the BH. Approximately 500 Mq of stars 
remain bound to the black hole. With these assumptions 
we do not account for stars that are originally unbound to 
the black hole. However, unbound stars do not contribute 
significantly in numbers deep within a full cusp. 

We run the simulations for 10^" yr, the approximate age 
of the clusters. Stars are removed from the simulation if they 
are ejected from the cluster, iJ > 0, if they reach a separation 
a > 10 pc, or if they are tidally disrupted by the central BIl(f|, 
^pcri ^ '"tid = Ki,{M,/mi,)^^^ . For each star that is tidally 
disrupted, we add its total mass to the BH. All of our A''- 
body simulations follow clusters with M, = 10* Mq and 
Vk = 5.6(T* ~ 105kms~^. For all of the reported simulations 
we set the time step criterion ry = 0.1 and the stars orbits 
were evaluated at minimum 80 times per orbit. We have 
checked the simulations with more precise parameters {rj — 
0.01 and 160 evaluations per orbit) and found similar results. 
We simulate 40 different cluster realizations for each Model 
and average over all the runs. A typical simulation takes up 
to one month on a single core of the Odyssey Cluster at 
Harvard University. 

2.2 A^-Body Results 

2.2.1 Cluster Evolution and Expansion 

Immediately after the recoil, the BH is ejected with 
~ 600 Mq of bound stars, in rough agreement with our 
initial estimates. In Figure [1] we plot the average number 
of stars bound to the recoiled BH as a function of time for 
Models I & II. After a relaxation timescale, ~ 10^ — 10^ yr, 
the star cluster begins to expand and lose stars as a power- 
law with Nci{t) (X t~^^^. Approximately 40% of the stars 
are ejected from the cluster and another 40% of the stars 
are tidally disrupted by the BH (see § I2.2.2p . After both 
Models begin to evolve, the ratio of the total mass of stars 
in each Model remains constant in time. 

In our simulations there are effectively only 3 param- 
eters that determine the cluster evolution: the BH mass, 
M, , the average stellar mass, m*, and the kick velocity 
Vk = 5.6(T*. All other scales in the simulation were deter- 
mined through the M. — a^, relatiorjf]. Since the required 
kick velocity to eject the BH scales with a*, and the to- 
tal number of stars in the cluster scales with the BH mass 
(cx M.), we can rescale the results of our simulations to find 
the number of stars in a recoiled cluster as a function of time 
{tr ocMf'''*) at t> lOVr, 



N,i{t) ^ 20 



M 







< m* > 



3/2 



f M. 



V 10* Mq 



13/8 



( * 



V 1010 yr 



-1/2 



(4) 



* We have also simulated clusters with a much smaller tidal dis- 
ruption radius in order to look at how this may affect the inner 
density profile. 

^ The tidal disruption radius depends on the mass ratio of the 
. phjIlH and star, however we found that the final number of stars 
in our simulations was insensitive to the chosen tidal disruption 
radius. 
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Figure 1. Evaporation of tlie star cluster in A'^-body simulations. 
The number of stars in the cluster bound to the BH is plotted as 
a function of time. The solid black line is for the cluster with the 
fiat mass function (Model II), the dotted red line is the cluster 
with equal- mass stars (Model I), and for comparison, the dashed 
blue line is Model I renormalized to have the same number of 
initial stars as Model II. After ^ 10^ yr the cluster decays as a 
power-law with time, TV cx t"^/^. Model I has approximately half 
the number of stars of Model II, because the stars are, on average, 
twice as massive. 



for a fixed Vk = 5.6(j*, and wliere we have normalized 
the scaling relation to match our A''-body simulations. The 
more massive the BH is, the fewer number of relaxation 
timescales the cluster will undergo over a fixed duration 
in time. Equations ((3]) & (|4]) are equivalent at t = 10^" yr 
when M. « 2 X 10® Mq. For BHs with masses above this 
value, the total number of stars in the cluster matches the 
initial number of stars at recoil. This has important im- 
plications for recoile d clusters near elli ptical galaxies and 
galaxy clusters (see iMerritt et all 120091 '). For star clusters 
with A4, < 2 X 10® M0 , such as those we expect around the 
Milky Way, the cluster will lose many of its initial stars and 
have evolved from its original state. From equation Q, and 
the power-law nature of Figure [1] for a fixed average stellar 
mass, we expect the final number of stars in the cluster will 
not be significantly larger even if the cluster had more stars 
initially; it would instead begin to decay earlier, following 
the same overall functional form N(t) oc t^^^'^. 

In Figure [2] we plot the normalized velocity spectrum 
of stars ejected from Models I and II as a function of time, 
looking at the first, second, and last third of stars ejected 
from the cluster. The typical velocity of an ejected star is 
usually a fraction of the velocity dispersion of the cluster, 
which decreases with time as the cluster expands. The to- 
tal velocity distribution of all the stars ejected from cluster 
appears almost log-normal, independent of the mass func- 
tion of stars, with a peak at ~ 6kms~^, and full-width half 
maximum an order of magnitude above the peak value. On 
average, nearly twice as many stars were ejected from Model 



II (~ 380 per cluster) compared to Model I (~ 220 per 
cluster). However, the total mass of stars ejected is com- 
parable. From the velocity spectrum, we see that the slow 
diffusion of stars to higher energies can not be the cause 
of the cluster evaporation. If this were the case, the ejec- 
tion spectrum would peak near zero velocity, as most stars 
become unbound just as they approach the escape veloc- 
ity of the cluster. In fact, the peak velocity of the ejected 
stars is near At; ~ t) ~ a, a reflection of large angle scat- 
tering (|Henonlll969l : iLin fc Tremaindll980l '). Indeed, this is 
confirmed by the pericenter distance of the stars before they 
are ejected, which is always much smaller than the half-mass 
radius of the cluster and the boundary of our simulations, 
'"peri << rfi « 10 pc. Nearly all of the ejected stars will 
remain bound to the Milky Way halo [v < 500kms^^). 
Only a few of the earliest stars ejected from the cluster have 
large enough veloc i ties to constitute a h ypervelocity star 
iBrown et al.ll2005l : IYu fc Tremainell2003l '). This, of course, 
is a small fraction of the number of hypervelocity stars that 
are exp ected to be produced during the inspiral of the BH 
binary (|Yu fc T remaine"2003l'). 

After approximately one relaxation timescale, the clus- 
ter begins to expand as it is heated by the innermost star 
in the cluster as well as by tidal disruptions. We find in our 
A'^-body simulations that radii that enclose a fixed number 
of stars scale as a power-law rjv cx t'^^^. The same relation 
is observed in our Fokker-Planck simulations (see § 13. ip , in 
previous simulations of bla ck holes in star clusters (see, e.g., 
lAmaro-Seoane et al.ll2003) . as well as in previous work ex- 
ploring the expansion o f a cluster wi t hout a black hol e post 
core-collapse (see, e.g.. lHenonlll96ll : lGoodmanlll984l . and 
citing articles.). The power-law index can be obtained by 
looking at the fiow of energy through the cluster, so long as 
the ene rgy is generated in a sufficie ntly small volume. Fol- 
lowing, [Oeierelal] iloil) (see also, lHenonlll96ll , ll965l ), we 
can analyse the fiow of energy 



E _ ^ r-N _ C 

E N VN fr(rjv) 



(5) 



at a radius that encloses a fixed number, N, of stars. Here, 
^, is independent of N, tn and E. Under the assump- 
tion that the rate of stars being ejected from the cluster 
is small {N/N « rjv/rjv), we can solve for rjv as a func- 
tion of t. For a cluster in the Keplerian potential of a black 
hole, tr(rjv) oc cr(rjv)^r^ oc r^^. Solving Eq. O then yields 
rjv oc t^''^. A similar relation can be found for the expan- 
sion ofa_£hister without a black hole post core-collapse (see, 
e.g.. lHenonlll96ll , [1965; Goodman 1984 . and citing articles.), 
which yields the same proportionality, however for a differ- 
ent reason {Uh oc N^^'^rf/^). As the ejection of stars from the 
cluster becomes important, however, we expect the cluster 
evolution to deviate from rjv oc t'^^^. 



2.2.2 Tidal Disruption of Stars & Resonant Relaxation 

A star will be tidally disrupted by the BH if it comes within a 
radius rtid ~ R*(M,/Mt)^'^''. Using this criterion, we remove 
stars that are disrupted by the black hole, and add their 
mass to the black hole. In Figure |31 we plot the average 
tidal disruption rate of stars in Models I fc II as a function 
of time. We find that after the first relaxation timescale, the 
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Figure 2. Velocity Spectrum of Ejected Stars. Plotted is the normalized velocity distribution of ejected stars from equal-mass clusters 
(solid black line, left panel) as well as from clusters with a fiat IMF (solid black line, right panel), Models I and II respectively. The red 
dotted, green dashed, and blue long-dashed lines represent the velocity distribution of the first (ejection < 10^ yr) , second (~ 10^ — 10* yr) , 
and last third (> 10* yr) of stars ejected owing to two-body scattering, respectively. As the cluster expands, the typical velocity of ejected 
stars decreases. The vast majority of stars are ejected with velocities less than the velocity dispersion of the Milky Way halo. The stars 
with the largest ejection speeds are ejected early in the cluster evolution, typically at < 10^ yr. 



time evolution of the disruption rate is well approximated 
as a power-law oc t~^^'^ . 

If a star has a small enough angular momentum such 
that its pericenter distance is less than the tidal disrup- 
tion radius, it will be disrupted in less than one orbital pe- 
riod. Therefore, the tidal disruption rate is determined by 
the flow of stars into the empty loss-con^. The stars can 
diffuse into the loss-cone through regular two-body relax- 
ation, large-angle scattering, or coherent e f fects such as res- 
onan t relaxation (|Rauch fc Tremainelll99d : lRauch fc Ingalld 
1 199^ ). We can hope to determine the primary mechanism be- 
hind loss-cone refilling using the time evolution of the system 
(i.e., Ntd oc ^ j^i^y 

The diffusion rate of stars into the empty loss-cone from 
two-body relaxation scales approximately as ~ N/tr- After 
the cluster begins to expand, the cluster should expand such 
that the relaxation timescale follows the clusters age, tr ~ t. 
Thus the rate of tidal disruptions from regular relaxation is 
oc t~^^'^, in agreement with our results. Large-angle scatter- 
ing will disrupt stars with a similar dependence on time, but 
at a rate reduced by ~ (lnA)~^. 

For resonant relaxation, the tidal disruption rate should 
be r; 'yN{< r)/trr{r), where trr{r) is the resonant relax- 
ation timescale, and 7 normalizes the rate and can be dctcr- 
mined us ing num erica l simula ti ons (Ranch fc Ingalla 1998 ; 
Hopman fc Alexa nder! l2006al : iKomossa fc MerritlT I2OO8I : 



Eilon et al.ii200S '). If we exclude general relativistic preces- 



^ The stars can also diffuse to higher specific energy, however the 
fractional change in energy required is usually much larger than 
the fractional change in angular momentum. 



sion, the resonant relaxation timescale is determined by the 
precession of stars from their own self-gravity, and scales 
roughly as trr ~ P{r)M,/m^, independent of the density 
profile of the stars. For a homologously expanding cluster 
around a black hole, the radius that encloses a fixed number 
of stars scales as oc t^^^ (see ^ I2.2.ip . Assuming that the or- 
bits are nearly Keplerian, P{r) oc r^^^'^ , the disruption rate 



from resonant relaxation will scale as N/ {t 



2/3\3/2 



N/t, 



the same as for regular relaxation. As can be seen in Fig- 
ure [S] the tidal disruption rate scales as N/t oc t~^^^ after 
t ~ 10^ yr, when the cluster begins to expand. From scaling 
arguments alone, we can not determine the relative contri- 
bution of tidally disrupted stars from resonant relaxation or 
regular relaxation. Resonant relaxation should be more im- 
portant on these scales since the total enclosed mass in stars 
is much less than the BH mass. However, artificial numerical 
precession can prevent resonant relaxation from occurring in 
numerical simulations. 

To determine the cause of the tidal disruption events 
we can analyze the tidal disruption rates dependence on m*. 
For resonant relaxation the disruption rate scales as oc m*, 
whereas for regular relaxation the rate scales as tr oc mj. 
Following Komossa fc Mcrritt (2008), we have performed 
smaller numerical simulations for recoiling star clusters with 
varying m* . These simulations had only ~ 200 of initial 
stars on an n oc density profile. This profile was cho- 
sen so that the regular relaxation timescale was shortest 
at largest radii. The tidal disruption rate as a function of 
time and mass is shown in Figure [l] Despite the stochastic 
variations given the small number of stars, we see that the 
tidal disruption rate scales approximately as oc at early 
times, and falls off as t~^^'^ after approximately one relax- 
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Figure 3. Tidal disruption rate of stars as a function of time. 
We plot the average tidal disruption rate of the star cluster as 
a function of time since it's central BH was kicked for Models 
I (single mass stars - dotted red line) & II (solid black line). 
The long dashed line is a power-law drawn for compari- 

son. For resonant relaxation the tidal disruption rate scales as 
N/trr, which is roughly constant at early times, and then scales as 
f-3/2 during the expansion of the cluster. Interestingly, although 
the intrinsic rate of tidal disruption events in a given cluster is 
small, the total rate of all recoiling clusters may contribute up 
to 10~® — lO"'^ yr~^, only one to two orders of magnitude lower 
than the disruption rate of black holes fixed in galactic nuclei. On 
the order < 10* yr, there is an initial burst of tidal disruptions 
partially because we started with a full loss cone in our simula- 
tions, however even if the loss cone is empty a burst may occur 
because the kick can put stars into orbits within the loss cone 
dStone fc Loebll2010h . 



ation timescale of the system (cx mj). We have confirmed 
that the disruptions are indeed caused by resonant effects 
by observing that the stars that are disrupted are preferen- 
tially from the inner most region of the cusp, and undergo 
angular momentum evolution on a timescale much shorter 
than the relaxation timescale. 

Because the timescale of large-angle scattering ejecting 
stars is also proportional to the relaxation timescale, both 
the tidal disruption rate and ejection rate have the same 
functional dependence on time. By the end of the simula- 
tions, approximately 40% of the stars are disrupted by the 
BH. T his is inconsistent with the results of lLin fc Tremain3 
who found that the BHs in globular clusters are more 
efficient at ejecting stars from density cusps than consuming 
them. However, here we are analyzing only a fraction of the 
entire density cusp. 

The tidal disruption of a star from an offset BH 
is an exciting possibilit y for detecting recoiled BHs 
l|Komossa fc MerrittI |2008| ). Previous work has so far fo- 
cused on the disruption of stars from clusters that were 
recently ejected by merger jKomossa fc MerrittI l2008( ). as 
it was thought that the tidal disruption rate would decline 



Figure 4. Tidal disruption rate as a function of stellar mass. We 
plot the tidal disruption rate of stellar clusters with = 10, 
1, and 0.1 Mq as green (dash-dotted), black (solid), and red 
(dashed) lines respectively. Initially, disruption scales approxi- 
mately as oc m^,. After approximately one relaxation timescale 
(cx m^) the cluster begins to expand, and the tidal disruption 
rate begins to decrease as oc t~^^'^. There is considerable scatter 
in these figures because we used only one simulation for each line. 



exponentially with time. This is in contrast with the power- 
law decline in tidal disruptions found here. Although the 
tidal disruption rate peaks early in the cluster evolution, 
most recoiled BHs are ejected from their host galaxy in the 
early Universe. Taking the time evolution of the clusters 
from equation Q, the tidal disruption rate of stars in a 
cluster today is approximately 



10" 



yr 



M. 



10'' M0 



13/8 



1010 yr 



-3/2 



(6) 



for each cluster with M, < 2 x 10^ Mq and t > tr. In 
Paper I, we found that there are perhaps tens of clusters 
with M. > 10^ M0, and hundreds with M. > lO" 
around each Milky Way like galaxy. Thus, per Milky Way 
galaxy, the tidal disruption rate of stars by rogue BHs 
is approximately 10"^ yr"i (10"'' yr^^) for M. > 10^ M© 
(> 10'* Mq). This is somewhat lower than the estimated 
tidal disruption rate of stars by BHs that reside in galax- 
ies ^ 10~^ — 10~^yr~i. Upcoming optical surveys, such 
as PTfQ, Pan-STARRfl and LSSTn are most sensitive to 



flares from BHs AI, 



10^ - 10^ Mq (Strubbe fc QuataertI 
|2009:), precisely the range of BHs we expect to have the 
highest present day tidal disruption rates. These events can 
be identified as off-nuclear flares, which have broad emis- 
sion lines with a mean redshift comparable to their nearby 



^ http : //www. astro . caltech.edu/ptf/ 
* http : //pan- Starrs . if a.hawaii . edu/ 



|http : //www . Isst ■ org/ | 
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galaxief^. Tidal disruption flares from low mass recoiled 
BHs may be a promising avenue for detecting these unique 
systems. 

2.2.3 Present State of Star Clusters 

The purpose of these long-term A'^-body simulations is to 
determine the present day distribution of stars around re- 
coiled BHs, with the goal of optimizing the search strategies 
for these clusters in § U In Figure [S] we plot the average 
projected number of stars enclosed within a radius r from 
the BH. In addition, we plot the results of our Fokker-Planck 
simulations from §[3] and Paper I, rescaled in the final num- 
ber of stars to match the A'^-body simulations. The stars 
in the A'^-body simulations are distributed with the same 
functional form and slope with a ~ 2.15 near the half mass 
radius, as we found in our initial simulations, despite the 
total number of stars being 1/10 of that found in Paper I. 
Regular relaxation, which was included in Paper I, appears 
to determine the shape and expansion of the cluster, whereas 
a combination of strong-encounters between stars and the 
tidal disruption of stars from resonant relaxation, neither 
of which were included in our Fokker-Planck simulations of 
Paper I, determine the final number of stars in the cluster. 
We conclude from this agreement that recoiled clusters of 
comparable mass stars have power-law density profiles with 
a < 2.15. 

We find only moderate evidence of mass-segregation in 
these clusters, even though the stars spanned a factor of ^ 10 
in mass, as shown in Figure [6] In our analysis of Model II, 
we binned the stars into shells around the BH, and found 
that the average mass of stars in each shell decreased as a 
shallow power-law of radius from 0.55 M© and 0.45 M©. Be- 
cause the massive stars are more luminous, this segregation 
will steepen the light density profile in the cluster. 

To check the robustness of our results on the underlying 
assumptions of our simulations, we have run additional sim- 
ulations that i) start with a shallower density profile, ii) in- 
clude stellar evolution, iii) include general relativistic effects 
to 2.5 post-Newtonian order, iv) have a significantly smaller 
loss-cone, or v) use a significantly different mass function. In 
cases (i) - (iv) there was no discernible effect on the stellar 
density distribution or final number of stars in the simula- 
tions. In case (v), where the average stellar mass was ten 
times larger, the cluster dissolved in less than 10^" yr. This 
simulation started with 190 stars in orbit around the black 
hole. At the end of the simulation, only two bound stars 
remained. 



3 FOKKER-PLANCK SIMULATIONS 

We found in our direct A''-body simulations that only ~ 10 % 
of stars initially bound to a recoiled lO** M0 BHs remained 
after 10^" yr. In contrast, our time-dependent Fokker-Planck 
simulations in Paper I, had little mass loss. Nevertheless, the 



Because the recoiling clusters are completely in the empty loss 
cone regime, the typical pericenter distance of the tidally dis- 
rupted star in a recoiled star cluster is always close to r^\^. This 
may result in a qualitatively different flare than associated with 
a central black hole. 



0.7 




0.001 



0.1 

Radius [pc] 



Figure 6. Mass segregation in recoiling clusters. Plotted is the 
average mass of stars as a function of radius when t = 6 X 10* yr. 
The average mass (with 100 stars per bin) slowly declines as a 
shallow power law out to the half-mass radius f» 0.2 pc. The 
lowest mass stars are preferentially scattered onto eccentric orbits 
with larger separation. 



Ai'-body simulations and our results from Paper I agree on 
the shape and slope of recoiled clusters. The A^'-body simula- 
tions show that strong encounters between stars as well as an 
enhanced tidal disruption rate drive the evaporation of the 
cluster on timescales much sh orter than standard, uncorre- 
lated perturbative encounters (|Henonlll969l : lLin fc Tremaind 
Il980l ). Here we reintroduce the time-dependent Fokker- 
Planck equation for st ars around a central m assive object, 
as originally derived bv lBahcall fc Wolj (Il976l ). with the ad- 
dition of two new sink terms in order to account for mass 
loss caused by stro ng encounters and resonant relaxation. 

Following (Ba hcall fc Wolj 119761 ) ■ we define the relax- 
ation timescale of the initial cluster to be 



3(271(7 



2\3/2 



327r2G2m2n* In A' 



(7) 



where cr* is the stellar velocity dispersion after the galaxies 
merge, m* is the average stellar mass, n* is the number 
density of stars at = GM./ctJ, and InA « ln(M,/Af*) is 
the standard Coulomb logarithm. In the dimensionless units 



of time, r = t/tr, and energy, x = —E/{mi,a-i, 
dependent Fokker-Planck equation reduces to 



the time- 



dg{x,T) 
dr 



5/2 d ^1 \ 



Rl,{x)~ Rrr{x)~- Rss{x), (8) 



where g{x,T) = [(2na'^)^^^n^^]f{E) is the dimensionless 
distribution function of the stars, Q{x) is the rate at which 
stars flow to higher energies, Ric{x) is the tidal disruption 
rate of stars that diffuse into the BH loss con e via reg- 
ular two-body relaxation (|Bahcall fc Wolj [l977l '). Rrv{x) 
is the tidal disruption rate of st ars that fall into the BH 
loss cone via resonant r elaxation llRauch fc Tremaind [199^ : 
iRauch fc Ingallj 1 19981 : iHopman fc Alexandej l2006al ). and 
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Figure 5. The total projected mass in stars within a distance r of the BH. Plotted is the cumulative mass projected within a circle of 
radius r for Model 1 (single mass: long-dashed red line) & Model II (flat IMF: solid black line) with A/, = 10* Mq. Fig. (a) shows the 
result of the fiducial A'-body simulation with the appropriate sized tidal disruption radius (see 12.111 . Overlayed, in the green dashed 
line, is the result of our Fokker-Planck simulations from Paper I chosen to have the same half-mass radius as our A^-body simulations and 
renormalized to have the same total number of stars. The overall shape of the curve matches the Af-body results remarkably well until 
there is, on average, less than one star enclosed. Also plotted (blue dotted line) is the result from § |3] including large angle scattering 
and resonant relaxation. The line was chosen at the time where the total mass of stars in the Fokker-Planck simulation matched the 
W^-body simulation of Model II. The Fokker-Planck simulations with resonant relaxation tend to clear out too many stars within the 
half-light radius, whereas the simulations without resonant relaxation have too many stars and must be rescaled in size in order to match 
the distribution of stars in the simulations. Fig. (b) shows a similar simulation of Model II with a much smaller tidal disruption radius 
^tid = 3iJs = 6GM, /c^ . The green dashed line in Fig. (b) doesn't include two-body scattering. The smaller r^i^ results in a density 
profile with power-law slope a = 2.15 throughout the entire cluster, where as, for a more realistic tidal disruption radius the density 
profile fiattens within the half-mass radius of the cluster. 



Raa{x) the rate that sta rs are ejected from the clust er ow- 
ing to strong encounters. iBahcall fc Woi3 (|l976l . [l977l ) found 
that 



Q{x) 



and 



d,[max(.,,)]-/M.(x)^-,(y).^^(^^ 



dy 



dx 



Rlcix) 



ln[J,(x)/JL 



(10) 



where Jc{x)/J-lc ~ {x/xtdY^"^ , and where Xtd ~ 
(M,/m*)~^'''^ri/i?* is the maximum specific energy of a star 
before tidal disruption. 

In systems with nearly Keplerian orbits, torques 
between the stars add coherently, and can efficiently 
randomize the angular momentum of the stars on a 
timescale much shorter than the regular (non- coherent) 
relajca tion time. This pro cess, known as resonant relax- 
ation l|Rauch fc Tremaine|[T996, ). can lead to an enhanced 
rate of tidal disruption s as the sta rs enter the loss cone 
Ranch fc Ingallj 1 19981 ). We follow iHopman fc Alexander! 



2006al ). who derived an approximate expression for the res- 



onant relaxation driven tidal disruption rate for Eq. [H] They 
found that the rate is approximately 



Rrr{x) ^ X- 



(11) 



where x is an unknown efficiency factor of order unity, and 
T„(x) is the resonant relaxation timescale. In our numeri- 
cal simulations we do not include general relativistic preces- 
siorF^. so the only form of precession that limits the reso- 
(9) nant relax;ation is caus ed by the enclosed mass of the s ystem. 
In these circumstances iHopman fc Alexanderl (|2006al ) found 
r„ ^ .0278a;3''2^ 

Although strong encounters are less important in cal- 
culating the flow of stars to higher and lower energies, in 
eq. (O, a single stro ng encounter can eject a star from the 
cluster. iHenonI (|l969 f) first calculated the escape rate of stars 
from isolated star clus ters for an arbitrary distribution of 
stars. iLin fc Tremainel (jl980i ) extended this work to stars 
around a central point mass. They found that strong en- 
counters are important in calculating the flux of stars out 
of the cusp, as confirmed in our A'^-body experiments in § [21 
By changing the limit s of integration in equation (35) of 
iLin fc Tremainel (Il980l '). we derive the rate that equal-mass 
stars are ejected from the cluster as a function of energy, 

,J,4,).3 V2 ^ P(^)d^ 1 (12) 

2 (x - xo)^ J (y + X - xo)^/^ InA 



In ii l2.2.3l we have included general relativistic precession and 
found no quantitative different in the tidal disruption rate or the 
distribution of stars in the cluster 
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where < is the negative specific energy required to be 
ejected. Note that in our dimensionless units, equation (|12p 
is suppressed by the Coulomb logarithm (In A)~^, compared 
to the rest of equation (|8]). 

We determine the time evolution of the cluster by nu- 
merically solving equation (|8]) with the boundary conditions 
g{x < 0) = exp(2;) and g{x > xtd) ~ until r = 2, at 
which point we set 17(2; < 0) = 0. We remove stars from 
the kick by scaling the distribution function of stars as 
g{x) g{x)z'^-^/{l + z^-^), where z = 2a;/(ufc/cr*)^. This 
yields an asymptotic density profile wit h n oc r"'* for r r^, 
as ex pected immediately after the kick (jKomossa fc Merritj 
l2008l ). We use a variety of xo = 0.01, 0.1, 0.25, 0.5, 1.0, 2.0, 
10.0, to explore the importance of xq in matching the num- 
ber of stars in §[11 as weU as x = 0.1, 0.5, 0.7, 0.8, 0.9, 1.0, 2.0 
to explore the uncertainty in the efficiency of resonant re- 
laxation. In all of our calculations we set InA = 10. 



3.1 Fokker-Planck Results 

In our calculations, we find reasonable agreement between 
the Af-body simulations and the numerical solution of the 
Fokker-Planck equation only when we include a new sink 
term, Rrr, which accounts for the tidal disruption of stars 
owing to resonant relaxation. We find the best agreement 
when we set the resonant relaxation parameter x ~ 0.8. 
When we exclude resonant relaxation, as we did in Paper I, 
we only get the proper functional form of the A'^-body so- 
lution to the density profile, but not the proper number of 
stars. Without resonant relaxation, the Fokker-Planck sim- 
ulation expand too rapidly at the half-mass radius. Indeed, 
when comparing the radii that enclose the inner 1, 10, or 100 
stars of the Fokker-Planck simulations to the A^-body simu- 
lations we find complete agreement that the cluster expands 
self-similarly as r oc t^^^ (see the discussion in ^ I2.2.ip . At 
the radius that encloses half of the cluster mass, we find the 
Fokker-Planck simulations still scale as t^^^ , and the A'-body 
simulations scales as t^^^. The t^^^ scaling is in agreement 
with our expectation that the relaxation timescale should 
be approximately the age of the cluster (Paper I). 

Resonant relaxation reconciles the two fundamental dif- 
ferences between our previous Fokker-Planck simulations 
and our current A^-body simulations. Resonant relaxation 
destroys stars at a faster rate than regular relaxation, re- 
ducing the number of stars in the cluster. As a consequence 
of the depleted stars, the outer cluster expands more slowly, 
matching the N-body simulations. In Figure [SI we plot the 
enclosed mass profile of the clusters from our TV-body simu- 
lations, along with our Fokker-Planck simulations including 
and excluding resonant relaxation (with x ~ 0.8). Resonant 
relaxation reduces the number of stars in the simulation to 
the correct value. Although including resonant relaxation 
in the Fokker-Planck simulations does give a similar mass 
density profile, there are too few stars in the inner most re- 
gion. Unfortunately, the Fokker-Planck simulations are one- 
dimensional, and can not take into account the anisotropy 
that must develop as the cluster expands, or the preference 
to deplete eccentric orbits near the black hole. 

The anisotropy of the cluster is often measured by the 
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Figure 7. Radial anisotropy of the expanding cluster. Plotted 
is the anisotropy parameter, /3 = 1 — uf/a^ as a function of 
radius for Model II when t = 6 X 10* yr. The cluster shows a 
large degree of anisotropy at nearly all radii. From the innermost 
stars outward, the stars in the cluster move from preferentially 
tangential orbits (/3 < 0) to radial orbits (/3 1). In our Fokker- 
Planck simulations we assume isotropy, i.e., /9 = 0, at all radii at 
all times. 




Time [yr] 

Figure 8. The fate of stars surrounding recoiled BHs. The num- 
ber of stars in the cluster bound to the BH is plotted as a func- 
tion of time for BHs with M, = lO'' Mq to 10* Mq from top to 
bottom. The evaporation of the clusters was calculated with the 
Fokker-Planck simulations including the tidal disruption of stars 
from resonant and regular relaxation (x = 0.8), as well as the loss 
of stars from strong encounters {xq = 0.1). For M« > 10® Mq, 
the cluster loses no more than ^ 60 % of its original mass. 
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parameter 



1 



(13) 



where ar is the radial velocity dispersion of the system, and 
at is the one-dimensional tangential velocity dispersion, such 
that crj = a"r + 2at . For an isotropic cluster, ar = at and 
P = 0. A cluster with stars on only radial orbits will have 
P — 1. We plot the anisotropy of the Ai'-body simulations as 
a function of radius in Figure [T] The cluster shows a large 
degree of anisotropy at nearly all radii. From the innermost 
stars outward, the stars in the cluster change from prefer- 
entially tangential orbits (/? < 0) to radial orbits (/? ~ 1). 
Indeed, the development of anisotropy is a natural conse- 
quence of the conservation of angular momentum in a Kep- 
lerian potential. As the cluster expands, the outermost stars 
must be on more radial orbits. The innermost stars, on the 
other hand, have their eccentric orbits depleted when they 
are tidally disrupted by the BH. In the our one dimensional 
Fokker-Planck simulations, we assume isotropy, which effec- 
tively relaxes the angular momentum of the stellar popula- 
tion on a timescale much shorter than the actual relaxation 
timescale. This causes the cluster to disrupt stars very close 
to the black hole, that would otherwise remain on circular 
orbits. This scenario can be tested wi th an appropriate two 
dimensional Fokker-Planc k code (e.g.. lCohnl 19801 : iTakahashil 
ll995l : lDrukier et al.|[l999h . extended to include the effects of 
resonant relaxation. Despite the anisotropy of the cluster, 
the density profile predicted by the Fokker-Planck simula- 
tions matches the A^-body simulations remarkably well, ex- 
cept in the innermost region of the cusp. 

The A-body simulations of § [2] are computationally 
challenging given the long timescale of the calculation and 
can not be easily extended to larger star clusters. Solving 
the Fokker-Planck equations, however, does not depend on 
Aci. Rather, it is calculated on a fixed grid in energy space, 
and takes only a short computational time to complete. We 
therefore use the Fokker-Planck code to compute the evolu- 
tion of more massive recoiled star clusters around M, = 10^, 
lO'', and 10^ M©, all with kick velocity scaled to their respec- 
tive velocity dispersion, Vk = 5.6cr*, and normalized by the 
M, — (T* relation. We set the free parameters x = 0.8 and 
xo = 0.1, which give comparable time evolution to the A- 
body simulations in § [21 In Figure [S] we plot the number of 
stars as a function of time for these recoiled clusters. We find 
that the evolution of star clusters around black holes with 
M. > 10® Mq, lose less than 60 % of their mass over 10^'' yr. 
For the largest black holes with M, > 10^ M©, the clusters 
lose negligible mass over the age of the universe. We can 
therefore expect that the most massive clusters represent 
the conditions of the stellar cusp when they were recoiled 
from their parent galaxy. 

We find that the total number of stars at the end of 
the simulations is only weakly dependent on the value we 
choose for xq. Indeed, over two orders of magnitude in xq, 
the final number of stars in the cluster changed by only 
10 — 20%. This is because for most of the range of xq values, 
the corresponding velocity for stars ejected with energy xq, 
was less than mean ejection velocity as seen in Figure (2] In 
the limit xq — > oo, equation (|12p goes to zero, and we recover 
the results of our simulations from Paper I. Unfortunately, 
we can not use these simulations to calibrate xq for BHs in 



galactic nuclei, where there is a reservoir of stars outside of 
ri. When xq ~ 0, equation (|12p is not accurate because it 
does not account for the fiux of stars to lower energy states 
from outside of the cusp or the return of stars tha t are still 
bound to the cusp of stars jLin fc Tremainelll980l ). 

The evolution of the star clusters is sensitive to the value 
of X- For sufficiently small, x ^ 0.1, clusters of stars around 
10* M0 BHs only lose ^ 60% of their mass over 10^" yr. 
Likewise, for x ^ 5, the cluster loses mass so rapidly, that it 
has < 1 star after 4 x 10^ yr. Such a scenario may represent a 
cluster with high concentration of stellar mass black holes. 



PHOTOMETRIC PROPERTIES OF OLD 
CLUSTERS 



SDSS DR7 (|Abazaiian et al.l 120091 ') has ~ 3.6 x 10* unique 
photometric objects that we would like to sort through in 
order to find ~ 100 candidate recoiled star clusters. Our goal 
is to develop the most general photometric model possible 
for recoiled star clusters in SDSS, while eliminating false 
positives as efficiently as possible. Generally, we expect the 
clusters to contain as few as ~ 20 stars for the smallest BHs 
M, = 10* M0 , up to ~ 10* stars for the most massive BH 
in the halo M, = 5 x 10^ M© . These clusters should have 
a power-law density profile with a ~ 2.15, but certainly 
1.75 < a < 2.25, which corresponds to a cusp of stars which 
flows away from the BH. Since the Milky Way has not had 
a recent major merger, we expect the clusters to be old. 
In this section, we develop such a model, focusing on the 
photometric properties of a stochastic cluster of old stars 
with a power-law density profile. Their spectra should indi- 
cate a large ct* for clusters around the most massive BHs, 
and large mass-to-light ratio at redshift z = 0, but we focus 
on photometric identification of candidates for spectroscopic 
follow-up. 

4.1 Cluster Models 

We generate model star clusters by randomly select- 
ing stars from a Kroup a initial mass function of stars 
dKroupa fc Weidneill2003l ) using two main modes of star for- 
mation. For Model A, we assume that the stars formed con- 
tinuously in time with a constant star formation rate until 
the BH was ejected after ~ 5 x 10^ yr. This is consistent 
with the e stimated star formation history of the MW galac - 
tic center (| Alexander fc Sternber^l 19991 : ICenzel et al.ir2003l '). 
We contrast this with Model B, where the stars formed si- 
multaneously with the BH merger, as galaxy mergers are 
often associated starbursts. In all instances we assume that 
the time between the merger of the galaxy and the merger 
of the BHs is negligible. 

The precise photometric properties of the stars depends 
on their metallicity and age. We consider three different 
metallicity histories of stars: (I) solar metallicity {Z — 0.02), 
(II) sub-solar metallicity (^ = 2 x 10~*), and (III) the esti- 
mated time evolution of metallicity of the galactic center. 

In all of our calculations we u se the online single star 
population synthesis code BaStQ l|Pietrinferni et al1l2004l : 



Available at |http : //astro ■ ensc-reimes ■ f r/basti/synth_pop/index ■ html] 
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ICordier et alJlioOTl ). and conv ert all observables t o the SDSS 
color system using Table 1 of IJester et alj (|2005l 'l. The typi- 
cal uncertainty in the colors from our conversion scheme is 
much smaller than the intrinsic variation in the colors of the 
star clusters. 

In Figure [9l we plot a color-color diagram of our Model 
Clusters with A^ci — 100, and teject ~ lO'^" yr after ejec- 
tion. Overall, the loci of star clusters follow the distribution 
of galactic stars identified by SDSS, however, they visually 
appear to be well separated from SDSS galaxies. In the fig- 
ure, it is clear that many clusters contain a star on the giant 
branch. Because such stars are more luminous than all of the 
other stars in the cluster, these systems are likely indistin- 
guishable from the population of late-type stars in the halo. 
Clusters with A'^ci = 10 — 1000 follow a similar distribution. 

Because the star clusters can not be distinguished from 
individual stars based on their colors alone, a successful pho- 
tometric search needs to exclude point sources, and use a 
magnitude system which doesn't depend on the exact light 
profile of the object. We therefore focus on the photometric 
properties of resolved objects, using the Petrosian m agni- 
tude system (jBlanton e t al. 2001; YasudaeFaDllooJ). This 
system is is defined by the Petrosian ratio. 



7^p(r) 



/p^^y dr'27rr7(rO/[7r(1.25^ - 0.8^)r^ 
/J'dr'27rr-'/(r')/(7rr2) 



(14) 



where I{r) is azimuthally averaged surface brightness profile 
in any particular band. The Petrosian radius, Vp, is defined 
by TZpifp) — 0.2 in the SDSS system. The total Petrosian 
flux (and magnitude) from the object is calculated as the 
total integrated light within 2rp, where rp is determined in 
the r-band alone. 

We use the azimuthally averaged cumulative light 
profile to distinguish candidates from galaxies and point 
sources. The SDSS catalog has the mean flux of light in 
annuli around the peak of the photometric object. We add 
the light in these annuli to recreate the total amount of 
light within radii n « 0.22,0.67,1.03,1.75,2.97,4.59,7.36 
arcsec. We use these bins to calculate the logarithmic 
slope of the cumulative light profile, Ti = dln/i/dlnr ~ 
ln(/i+i/Ji)/ In (ri+i/ri). For a completely resolved star clus- 
ter with power-law density profile n oc r"" , we expect the 
light profile to approach 3 — a ~ 0.85 where a ~ 2.15. Out- 
side the half-light radius of the cluster, however, the slope of 
the profile changes significantly from a single power-law. To 
model this we use the distribution of stars from our TV-body 
simulations and generate mock light density profiles with a 
variety of PSFs and angular sizes. To facilitate the search in 
§ [3 we have split our mock profiles into bins of varying Vp. 
The best fit parameters are detailed in Table [T] 



5 PHOTOMETRIC SEARCH 

SDSS has imaged approximately one quarter of the sky to 
a limiting magnitude r « 22.2. As the largest database of 
photometric and spectroscopic objects in the sky , it presents 
a prime opportunity to search for recoiled clusters. In § 14.11 
we developed a simple photometric model of recoiled star 
clusters in Milky Way like galaxies. Here we use this model 
to systematically search for photometric candidates in SDSS 
DR7. 



Table 1. Best fit SDSS light profiles to find mock clusters. Be- 
cause the cumulative light profiles are not corrected for seeing, we 
search for candidates by focusing on a range of properties that de- 
pend on the observed Petrosian Radius, rp. See S l4.1l for detailed 
definitions of the parameters. ''Seeing better than 1.2". ''Seeing 
between 1.2 and 1.7". 



r4 



re 



(arcsec) 



2.0-3.0=" 
2.0-3.2^ 
3.0-4.5 
4.0-6.0 



0.78 - 0.! 
> -85 



> 0.05 
0.25 - 0.5 
0.35 - 0.70 



0.10-0.25 

0.28 - 0.45 0.1 - 0.2 



We use the SDSS DR7 CAs,]OBf0 photometry database 
to select objects by size, shape, color, and azimuthally av- 
eraged light profile. We limit our search to resolved objects 
with a Petrosian radius Vp > 2" in the g band. 

Our color selection criteria is illustrated by the trape- 
zoids in Figure [9] We use the Petrosian magnitud e sys- 
tem color corrected for extinction (ISchlegel et al.|[l998l ) . Our 
criteria focus on an old population of metal-poor to so- 
lar metallicity stars but excludes clusters with a star on 
the red giant branch. We choose our candidate clusters 
out of the parallelogram defined by 1.25 < u — g < 1.75, 
0.5{u - g) - .225 < g -r < 0.5(« - g) - 0.075. Addition- 
ally, we require that the shape of the candidates be circular 
by selecting for objects with ratio of semi-minor to semi- 
minor axes greater than 0.7. These criteria select ~ 70, 000 
resolved photometric objects as candidate recoiled clusters, 
with photometric properties of stars. We limit our sample 
further by using the azimuthally averaged cumulative light 
profile (Ti) of the remaining objects as detailed in Table [1] 
In addition we search for simpler model clusters with r4 and 
Fs between 0.6 and 0.9. 

Using these criteria, we are left with ^ 1 , 000 candi- 
dates, which we visually inspect to remove obvious interlop- 
ers. These tend to be individual or binary stars in crowded 
fields, face-on disk galaxies, and cuspy elliptical galaxies. 
We are left with ~ 100 objects, which we list in Tables [2] 
& [31 Thumbnails of a selection of candidates is shown in 
Figure [TOl The number of candidates that remain is likely a 
reflection of our ability to visually inspect the candidates. It 
is impossible to visually inspect the 70, 000 objects selected 
through color alone, and we had to use some model depen- 
dent choices for the light profile to obtain a more reasonable 
number of photometric candidates (~ 1, 000). 



6 SPECTROSCOPIC SEARCH 

The selection of spectroscopic objects in SDSS DR7 is not 
ideal for serendipitously locating the star cluster around a 
recoiled BH. Indeed, many of the main science objectives for 
spectroscopic targets specifically exclude objects with the 
photometric properties similar to the recoiled clusters we 
search for. Nevertheless, resolved recoiled clusters would be 
identified in SDSS photometrically as a galaxy because of its 



Located at |http : //casjobs . sdss ■ org/CasJobs/"| 
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Figure 9. Color-color diagrams of model recoiled clusters. Plotted is the distribution of recoiled star clusters with Nd = 100 at 
Reject ~ 10^" yr with Z = .02 (Black Crosses), Z = .0004 (Red Crosses), and Varying Z (see text; Blue Crosses), along with a random 
selection of galaxies (magenta points) and stars (green points) from SDSS. The trapezoids correspond to our color-color selection (see 
§0. 




Figure 10. Thumbnails of a diverse selection of candidates. From left to right and top to bottom these are SDSS 
J114607.52-I-135233.1, SDSS J130154.22-031323.3, SDSS J052222.65-013302.9, SDSS J084034.69+162319.5, SDSS J084822. 47+355630.4, 
SDSS J093815.82-I-231234.8, SDSS J121414.73-f 161215.4, and SDSS J123544.93-I-193016.9. The scale is the same for all images, with the 
photometric object located at the center. 



extended size, but unlike galaxies, would have an extremely 
low redshift, z < 10'^. 

We use SDSS spectroscopy to select candidate clus- 
ters by their identified redshift with -0.002 < z < 0.002, 
corresponding to radial velocities with magnitude less than 
600kms~^. To exclude single and unresolved binary stars, 
we restrict our results to objects with Petrosian radii, rp > 
3.0 " in both the r and g bands. We then remove blended ob- 



jects, which are mostly galaxies with foreground stars. This 
results in approximately 270 objects identified by SDSS pho- 
tometry as galaxies, and 18 identified as stars. We have vi- 
sually inspected all 290 objects. All of the remaining objects 
in the sample have two main features: (i) a star like object 
on a diffuse source or (ii) featureless and diffuse source. In 
case (i), the spectral identification of the source is always 
stellar. The majority of objects that fall into category (ii) 
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Table 2. Candidate recoil clusters based on the selection criteria 
described in §[5] and Tab.[T] 
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were incorrectly identified with z = 0. The remaining cases 
were not spherical in shape or clumpy, and therefore not 
classified as recoiled cluster candidates. 



7 LITERATURE SEARCH 

We have also explored the literature for galaxies that were 
spect roscopically ide ntified as stars owing to their low red- 
shift. ISargentI ()l970l ) took low-dispersion spectra of 141 ob- 
jects selected from Zwicky's catalogs of compact galaxies 



Table 2. continued. 



Object 






; 


) 




7 


"p 


SDSS 


J120648, 


21+450646.7 


21 


.06 


20, 


.08 


3, 


,16 


SDSS 


J121700, 


,34+353542.1 


20, 


.38 


19, 


.74 


3, 


,22 


SDSS 


J123544, 


,93+193016.9 


20, 


.55 


19, 


.87 


2, 


.84 


SDSS 


J123614, 


,93+013708.7 


20, 


.27 


19, 


.81 


3, 


.03 


SDSS 


J125011, 


62-021800.1 


21 


.16 


20, 


.51 


2 


.65 


SDSS 


J125734. 


,92+253916.5 


20, 


.49 


19, 


.92 


2, 


.78 


SDSS 


J125915, 


28+070342.6 


20 


.01 


19 


.52 


3. 


.30 


SDSS 


J125958, 


,63-002508.4 


19, 


.75 


18, 


.71 


3, 


.26 


SDSS 


J130109, 


32+462607.1 


21 


.28 


20 


.15 


2, 


.76 


SDSS 


J130153, 


77+504842.2 


20, 


.38 


19, 


.64 


2, 


,85 


SDSS 


J130154, 


22-031323.3 


19, 


.67 


18, 


.97 


3, 


.32 


SDSS 


J132249, 


56+084115.6 


19, 


.85 


19, 


.22 


3, 


.27 


SDSS 


J132610, 


71+535511.7 


20, 


.79 


20, 


.26 


3, 


,24 


SDSS 


J133057, 


,40+184836.7 


20, 


.00 


19, 


.28 


2, 


,66 


SDSS 


J133212, 


59+353159.7 


19, 


.56 


18, 


.94 


2, 


.93 


SDSS 


J134127, 


13+081550.6 


21 


.56 


20, 


.52 


2, 


.06 


SDSS 


J134459, 


,96+030428.6 


19, 


.76 


19, 


.55 


3, 


.13 


SDSS 


J134737, 


53+203427.0 


20, 


.19 


19, 


.78 


2, 


.94 


SDSS 


J134852, 


16+245743.4 


19 


.72 


19 


.16 


3. 


.77 


SDSS 


J135040, 


,46+103538.7 


20, 


.75 


19, 


.95 


2, 


.84 


SDSS 


J135241, 


53+121430.8 


20, 


.51 


19, 


.66 


3, 


,06 


SDSS 


J135544, 


,47-065531.4 


20, 


.29 


19, 


.54 


2 


.90 


SDSS 


J140113, 


91+060627.7 


20 


.13 


19 


.71 


2, 


.85 


SDSS 


J141327, 


28+282847.1 


19, 


.58 


18, 


.79 


3, 


.05 


SDSS 


J141418, 


26+454312.8 


18 


.94 


18 


.57 


4, 


,37 


SDSS 


J142920, 


56+261616.5 


20, 


.13 


19, 


.39 


2, 


.80 


SDSS 


J142935, 


43+073722.6 


19 


.59 


18 


.87 


3. 


.18 


SDSS 


J145030, 


79+380441.6 


19 


.62 


19 


.06 


2, 


.91 


SDSS 


J145145, 


53+103402.0 


19 


.94 


19 


.32 


3, 


.45 


SDSS 


J145150, 


,02+352929.8 


20, 


.03 


19, 


.53 


2, 


,78 


SDSS 


J145345, 


50+080808.7 


20, 


.37 


19, 


.58 


3, 


,03 


SDSS 


J150113, 


,32+051304.1 


19, 


.92 


19, 


.37 


2, 


.91 


SDSS 


J150459, 


72+081819.7 


20, 


.50 


19, 


.89 


2, 


.78 


SDSS 


J151934. 


,33+134102.8 


20, 


.51 


20, 


.12 


3, 


.33 


SDSS 


J152006, 


,92+085031.0 


20, 


.45 


19, 


.68 


2. 


.78 


SDSS 


J152249, 


28+473700.1 


20, 


.31 


19, 


.82 


3, 


.12 


SDSS 


J152646, 


,00+210607.1 


21 


.21 


19, 


.97 


3, 


.33 


SDSS 


J153145, 


,65+150057.6 


20, 


.12 


19, 


.54 


3, 


,11 


SDSS 


J155333, 


,07+423146.0 


20, 


.25 


19, 


.53 


3, 


,07 


SDSS 


J155442, 


55+055111.1 


19, 


.09 


18, 


.60 


3, 


.23 


SDSS 


J160630, 


13+351046.2 


19, 


.48 


19, 


.06 


3, 


.04 


SDSS 


J160702, 


,48+110353.3 


20, 


.27 


19, 


.78 


3, 


.89 


SDSS 


J162536, 


57+563531.9 


20, 


.99 


19, 


.96 


3, 


.53 


SDSS 


J163339, 


,07+132635.6 


20, 


.23 


19, 


.56 


3, 


.46 


SDSS 


J163659, 


,29+235816.2 


20, 


.05 


19, 


.01 


2, 


.93 


SDSS 


J170525, 


,39+235241.5 


19, 


.17 


18, 


.37 


4. 


.08 


SDSS 


J172243, 


,89+080447.8 


20, 


.44 


19, 


.97 


2. 


.91 


SDSS 


J210803, 


13-001350.4 


20, 


.35 


19, 


.52 


3, 


.15 


SDSS 


J213035, 


,54-070545.7 


20, 


.76 


20, 


.05 


3 


.19 


SDSS 


J215424. 


,98+002023.4 


20, 


.01 


19, 


.14 


3, 


,23 


SDSS 


J233106. 


,11+075810.9 


20, 


.89 


19, 


.99 


2, 


.97 



(available in lZwickv fc Zwickvl[l97lh . and found that 14 ob- 
jects had near-zero redshifts and identified the galaxies as 
having a foreground star. We have reanalyzed newer digital 
images and spectra of these objects. In some cases the fore- 
ground star has moved and new spectra show the objects 
are extragalactic. However, most were visually identified as 
galaxies because of their disk like shape. Only one object, 
IV Zw 26, could not be excluded as a candidate owing to 
insufficient resolution in any survey. 



© 0000 RAS, MNRAS 000, 000-000 



14 O'Leary & Loeb 



Table 3. Candidate recoil clusters based on an asymptotic cu- 
mulative light profile with slope between 0.6 and 0.9. 
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8 SUMMARY AND CONCLUSIONS 

We have followed the long term evolution of star clusters 
around recoiled BHs using long term A^-body simulations, 
with a one-to-one correspondence between the stars and A^- 
body particles. We have found that for M, = 10'' Mq, ~ 40% 
of the stars are ejected from the star cluster, and ~ 40% 
of the stars are tidally disrupted by the central BH within 
10^° yr. We have scaled these results to BHs with masses 
M. < 2 X 10'' Mq, finding that iV,i = 840(M./10^ Mq)^^/** 
stars remain around the BH today for a typical recoiled BH. 
For more massive BHs, the cluster should eject or disrupt 
few stars over ~ 10^" yr. Although a single BH has a small 
tidal disruption rate, we have found that the total rate for 
all clusters in Milky Way like galaxies is ~ 10~^ yr~^, which 
is only a factor of ~ 10 lower than expected in the galac- 
tic center. We have extended our one-dimensional Fokker- 
Planck treatment in Paper I to include resonant relaxation 
and large-angle scattering to account for the dominate mass 
loss mechanisms of the cluster. Using this treatment, we 
were able to get satisfactory agreement between the A^-body 
simulations and Fokker-Planck simulations. Some discrep- 
ancy remained, which we attribute to the large amount of 
anisotropy in realistic clusters. 

We used our A''-body simulations to generate random 
realizations of star clusters today, which guided our search 
for star clusters around recoiled BHs in SDSS. In our pho- 
tometric search through SDSS data, we assumed that the 
star clusters have a power-law density profile and that they 
have colors comparable to an old population of stars. We 
used these criteria to find --^ 70, 000 candidates of which 
only ~ 1000 had a light density profile out to 4" consistent 
with a recoiled star cluster with power-law density slope. 
We visually inspected all candidates, and found that many 
were the bulges of nearby face-on spirals. The remaining 100 
candidates were faint, and difficult to distinguish from dis- 
tant galaxies. Follow-up spectroscopy is necessary to identify 
their nature. If any of them are a star cluster around a re- 
coiled BH, it would show unusually high velocity dispersion 
cr* at z ~ 0, with a spectrum of a population of old stars. 



We also searched the spectroscopic database of SDSS 
for resolved objects with a low redshift/blueshift consistent 
with the Local Group. The vast majority of these candidates 
were galaxies with bright foreground stars. We found no 
candidates that appeared to be a recoiled cluster. 

The criteria we used to search for candidate clusters 
required the cluster to have a well defined power-law den- 
sity profile, as found in our numerical simulations. Because 
we did not include compact remnants in our simulations 
(which would have the correct density profile but not light 
profile) we can not be confident that we properly included all 
star clusters in our search. Unfortunately, most research on 
stellar remnants around supermassive black holes focus on 
mass s egregation aro und relax ed stell a r cusps llFreitag et al] 

O'Leary et al.ri200a 



20061: iHopman fc A lexander l2006bl: 



Alexander fc Hopmanii200& : iKeshet et al.|[2009l ). These sim- 



ulations find that compact remnants play an important role 
in the dynamics for radii r <rk- We have not included these 
objects because it is difficult to quantify how many com- 
pact remnants would be in this region immediately before 
the binary merges. Indeed, the segregation of the compact 
remnants c an only occur after the low-mass stars populate 
this region (|Merrittll2009l ). In future work we hope to include 
these compact remnants, and extend our Fokker-Planck code 
to inclu de large angle scatt ering between stars of multiple 
masses (jO'Learv et al.ll2009l '). 

An alternative search strategy, which we did not explore 
here, is to cross-check less stringent criteria with alternative 
databases, such as the ROSAT x-ray survey. 
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